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FUNDAMENTAL ANALYSIS OF THE FAILURE OF 
POLYMER- BASED FIBER REINFORCED COMPOSITES 


INTRODUCTION 


There are many different ways in which a structure made of a fiber 
reinforced composite material can become unable to adequately perform its 
primary function. In each such instance failure is considered to have oc- 
curred. The possible failure modes therefore encompass a wide range of pos- 
sibilities from simple loss of structural rigidity due to gross inelastic de- 
formation (e.g., yielding), through a reduction in load-carrying capacity 
due to localized damage or separation (e.g., interply delamination), to the 
complete loss of strength due to large-scale crack growth and fracture. Each 
of these failure modes can be gradual or rapid depending on the nature of the 
applied loads, the material properties, the geometry of the structure, and the 
presence of cracks or flaws. For polymer-based composites, the loading rate, 
the temperature, and previous load history can also play prominent roles. 

In the work described in this report, emphasis is placed upon 
fracture and, therefore, the work will be primarily concerned with the "strength” 
of fiber composites containing known flaws. The term strength is convention- 
ally taken to mean the load level at which failure occurs in a standard test 
specimen. Clearly, the strength will be a function of many different para- 
meters arising in the test program and may or may not be directly applicable 
to engineering design situations. The primary purpose of the research des- 
cribed in this report is to provide a bridge between standard laboratory test 
procedures and actual engineering applications of fiber composites that will 
allow accurate reliable estimates of the failure loads for aircraft and other 
engineering structures. Such a capability does not presently exist. 

The design tool that should be developed for the safe and efficient 
utilization of composite materials is a predictive capability for fracture 
that can take account of the applied loading, the geometry of the structure, 
and the environmental effects in terms of readily measurable properties of the 
composite’s constituents and of its microstructural design. The primary 
benefits accruing from such an analytical capability are twofold. First, 
for a given structure, material, and load, the critical flaw size and type in 
a structural component can be estimated for comparison with actual flaws ob- 
served in an inspection program. Second, guidelines can be provided tor tne 
designer to tailor more fracture resistant "tougher" materials for specific 
engineering applications. 
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PROGRAM PLAN AND SUMMARY OF PROGRESS 


The work described in this report represents the first year of effort 
in what is planned as a multiyear program. The specific objective of the 
first year’s work was to develop the basic mathematical procedures required 
for the analysis model. To implement this work, attention was focussed on 
unidirectional composites with linear elastic-brittle material behavior. In 
subsequent work, the model will be extended to treat angle ply laminates and 
will include further refinements (e.g., inelastic constituent behavior) re- 
quired in order to treat actual engineering problems. 

The primary objective of the work will be achieved only if the 
mathematical model developed is capable of delineating the role of the various 
micromechanical failure processes that dictate the ultimate failure point of 
fiber reinforced composites. The research described in this report seeks this 
end by merging a micromechanical failure analysis with a macromechanical frac- 
ture mechanics approach. This approach treats the material as heterogeneous 
and anisotropic where microstructural effects predominate and as homogeneous 
and anisotropic where it is permissible and practical to do so. In this way, 
direct consideration can be given to 

• The external size and shape of the structure and the 
laminate stacking sequence 

• The applied loads acting on the structure, both 
mechanical and thermal, and environmental effects 

• The size, shape, and orientation of a flaw in the 
laminate. 

In particular, the manner in which these parameters influence the sequence 
of microstructural failure events whereby a flaw extends stably under an 
increasing load up to the point of catastrophic fracture will be determined. 

It should be understood that the conventional approach to fracture analysis, 
linear elastic fracture mechanics, is not capable of coping with this degree 
of complexity, thus necessitating the more general development being pursued 
here . 
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Linear elastic fracture mechanics (LEFM) is a predictive technique 
applicable to structural components containing crack-like flaws when the 
material used fits certain key assumptions used in the LEFM theory. Specifi- 
cally, the material must behave very nearly as would a completely linear 
elastic perfectly homogeneous ideal material. The LEFM approach has been 
successful when properly applied (e.g., to high strength steel), but consider- 
ably less successful in applications where the basic assumptions are not well 
satisfied. Fiber reinforced composite materials are an outstanding example 
of the latter case. In particular, in fiber composites with a flaw, a non- 
linear "damage" zone is generally produced at the flaw. This zone, generally 
growing in a stable manner under an increasing applied load, has a profound 
effect on the eventual point of complete fracture. This is shown, for example, 
in the work of Mandell, et al, [l]* as well as by many other investigators<, 

The damage zones in a fiber composite are the result of a large 
number of discrete failure processes, e.g., fiber breaking, matrix yielding, 
etc. , that occur in the highly stressed region ahead of a crack-like flaw. 

These individual micromechanical events do not conform to the basic assumptions 
of LEFM. Consequently, it is not surprising that the agregate of such pro- 
cesses cannot be treated by LEFM. What is therefore needed is a generalized 
fracture mechanics treatment which explicitly recognizes the fundamental na- 
ture of damage growth in composite materials. Specifically, a proper analysis 
of fracture in fiber composites must be cognizant of two key features: the 

generally anisotropic constitutive behavior of the material, and the hetero- 
geneous nature of the fracture process. 

The approach developed in this report can be likened conceptually 
to boundary layer theory and, in application, to the well-established singu- 
lar perturbation and matched asymptotic expansion techniques of fluid 
mechanics. That is, the problem of a composite material containing a flaw 
is divided into distinct "inner" and "outer" regions. In each of these regions, 
the material is modeled in different ways. The inner region, which contains 
the crack tip, is considered on the microscopic level and treats the material 
as being heterogeneous. The outer region surrounds the crack tip region. It 
is treated as a homogeneous ortho tropic continuum. For simplicity, the inner 
crack tip region will be referred to in this report as the LHR (for local hetero- 
geneous region) while the outer region will be simply called the continuum. 


* References are given on page ^^5. 
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The LHR consists of elements representing the matrix, the fibers, 
and the fiber-matrix interface zones. The constitutive relations of these 
elements, up to and including their rupture points, are presumed to be known 
from experiments. Any element of a fiber composite ruptures when an intrinsic 
critical energy dissipation rate can be provided at some point of that element. 
These critical values are assumed to be independent of the local stress field 
environment at the point of incipient rupture. This permits data from fracture 
tests on isolated fibers, matrix material, and, possibly, unflawed composites 
(to obtain interface strengths) to be directly inserted into the model- Ma- 
terial properties used in the analysis work and, ultimately, critical tests 
of the predictions of the model will be obtained from a concurrent experi- 
mental program being carried out under NASA sponsorship at Virginia Polytechnic 
Institute and State University. 

Progress to date has permitted computations to be performed for 
unidirectional composites with elast ic-perfectly brittle constituent 
behavior. The mechanical properties have been those of graphite epoxy. The 
rupture properties have also been arbitrarily varied to test the capability 
of the model to reflect real fracture modes in fiber composites. It has 
been shown that fiber breakage, crack bridging, matrix-fiber debonding, and 
axial splitting can all occur during a period of (gradually) increasing load 
prior to catastrophic fracture. In this way, the sequential and interrelated 
manner in which each individual local rupture event occurs during damage 
growth preceding fracture in fiber composites is revealed by the computations. 
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DESCRIPTION OF THE ANALYSIS TECHNIQUE 

In the preceding section of this report, the objectives of the work 
and an outline of the general approach were given. The focal point for this 
description was the two-dimensional local heterogeneous region (LHR) surrounding 
the crack tip. A typical LHR model is shown in Figure 1. Depicted is a uni- 
directional fiber composite containing three distinct components: the fibers, 

the matrix, and the fiber-matrix interface zones. In this section of the report, 
the discrete elements that comprise the LHR and the manner in which they are 
assembled to give a quantitative predictive capability for composite fracture 
will be described in detail. 


LHR Boundary Conditions 


In the work performed so far, the interaction between the LHR and the 
continuum, as is appropriate for preliminary stages of the work, has been taken 
in the simplest possible manner. This is by specifying the displacements on the 
boundary of the LHR in accord with the "rigid boundary conditions" approach. Use 
of this scheme is tantamount to assuming that the LHR is large enough that the 
nonlinear inhomogeneity of the crack tip region does not affect the periphery 
of this region. In other words, the LHR and the continuum are assumed to be 
uncoupled. Consequently, the displacements at the LHR boundary are exactly the 
same as if the entire cracked body is an elastic continuum. The required rela- 
tions can therefore be obtained from the work of Sih and Liebowitz [2]. Their 
approach is briefly summarized, as follows. 

Consider a plate with its major dimensions lying in the xy plane. 

For either plane stress or plane strain deformation, the inplane strains 

c , and Y will depend only on the inplane stresses a , a , and t 
y ' xy X y xy 

Hence, the appropriate constitutive relations for a rectilinearly anisotropic 
body in a state of plane deformation are given by 


e = o + 
y 12 X 

' xy 16 X 


+ a 


16 


+ a 


26 


+ a 


66 


T 

xy 

T 

xy 

T 

xy 


( 1 ) 
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FIGURE 1. THE LHK CONCEPT 
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Sih and Liebowitz show that the solution to the governing differential equation 
of two-dimensional anisotropic elasticity theory is associated with the roots 
of the characteristic equation 


‘^11 ^ ^66^ - 2 a2g 4 + a^2 = 0 . ^2) 

It can be shown that the roots of Equation (2) are either complex or are pure 

imaginary. Hence, these can be labeled 4 ^, ^25 -6^, and This suggests the 

introduction of the complex variables z_ = x + ^ y and z = x + ^^y. The plane 

1 1 2 ^ 

problem of an anisotropic body is thereby reduced to the determination of two 
complex potential functions of a complex variable cj)(z^) and ^(^ 2 ^ that satisfy 
the prescribed boundary conditions of the problem. 

As further shown in reference 2, the displacement field is expressed 
in terms of the potential functions by the relations 


u = 2 Re {p^ (j) (z^) + P 2 '|j (.Z^) } 


(3) 


V = 2 Re {q^ (() (z^) + q2 >1^ (z^) } 

where u and v are the x and y displacement components and 
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12 


- a 


16 


- “ii ^ 


+ a 


12 


- a 


16 "^2 


^^12 "^1 ^22 “^26 


2,2 ^“^12 ^2 ^^22 ^26 ^2 ^ 


Omitting the details, potential functions for a cracked body infinite in extent 
have been determined for insertion into Equations (3). In particular, a solution 
for remote loading consisting of a uniform tensile loading a acting in the di- 
rection normal to the crack plane and a shear loading t parallel to the crack 



8 


plane for a crack of length 2a can be obtained. For a polar coordinate system 
with origin at the crack tip as shown in Figure 2, the displacements near the 
crack tip are given by 




1/2 

^1 p 2 (cos (() + ^2 sin (j)) + 


^2 ?! (cos (j) + sin ([>) 


1/2 


K 


1/2 




> 1/2 


(cos (j) + ^2 sin (p ) + 


- (cos (}i + sin <j> ) ' 


(4) 


and 



(cos (p -h 


sin (p ) 


1/2 ^ 


” ^2 qj^ (cos (p sin <p ) 


1/2 


+ 


-II (^r 


1/2 

q 2 (cos (p 6^ sin (p ) + 


- (cos (j) + sin (p ) 


1/2 


I ■ 


(5) 
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FIGURE 2. DISPLACEMENT COMPONENTS FOR A POINT IN THE 
VICINITY OF THE CRACK TIP 



10 


where 


= a / TT a 


and 


( 6 ) 


= T /tt a 

are the Mode I and Mode II stress-intensity factors for the problem. 

Equations (4) and (5) can be used to determine displacements on the 
LHR boundary via the rigid boundary-condition approach. This means that the 
applied stresses acting on the body are transmitted through the continuum region 
to the crack-tip region and are "sensed" at the crack tip in terms of the stress- 
intensity factors and Kjj* An independent specification of the load and 
crack length is therefore unnecessary. Hence, although derived for an infinite 
medium, the approach can be used for bodies with finite boundaries by simply in- 
serting the appropriate stress-intensity factors. In doing this, it must be 
tacitly assumed that the LHR is (1) large enough relative to the micros true tural 
dimensions of the composite that the boundary displacements are closely given 
by continuum theory and (2) small enough relative to the crack length and di- 
mensions of the body that the singular behavior of the continuum solution at the 
crack tip dominates. 

While appropriate for preliminary work, it is anticipated that the 
rigid boundary condition approach — even with periodic updating to reflect the 
progress of the crack through the LHR — may prove to be too restrictive. There- 
for, in subsequent work, a flexible boundary-condition approach which extends 
an approach developed in work previously carried out at Battelle [3] will be 
used. This is described in more detail in the Recommended Further Research 
section of this report. It might be noted that with flexible boundary conditions, 
crack length and load must be specified individually. 
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The LHR Element 


As shown in Figure 1, the LHR for a fiber reinforced composite is 
considered to be made up of three different types of constituents: fiber, 

matrix, and fiber-matrix interface. Each of the individual constituents in 
the LHR must be capable of rupturing to allow the body to exhibit the changes 
in strength that correspond to various levels and orientations of local break- 
age. A good deal of success has been obtained by Kanninen [4,5] using "spring- 
like” elements to model a local fracture phenomena. This fact, taken together 
with the increased computational simplicity of this approach, has led to the 
adoption of the basic element shown in Figure 3 for this work. 

As shown in Figure 3, each element has eight degrees of freedom 
and is connected in the LHR at its four corner node points. Extensional stiff- 
ness is provided by four extensional connectors. These connectors resemble 
simple extensional springs, but also have a lateral contraction or Poisson effect. 
A material made up of a set of these connectors behaves in uniform extension 
exactly like a homogeneous orthotropic material. Similarly, shear stiffness 
is incorporated in the element through a rotary spring at each node point to 
give the proper response to shear loadings. The values of the spring con- 
stants are functions of the material’s elastic constants and of the element 
size and shape. 

In addition to giving the proper response to loads, the LHR elements 
fracture according to the "codes" shown in Figure 4. In Figure 4, Fracture Code 
1 represents an increment of crack extension in the plane lying midway between 
Nodes 1 and 4 that extends from the left side of the element to the center of 
the element. Fracture Code *2 represents the continuation of the crack to the 
right side of the element. Fracture Codes 3 and 4 represent similar increments 
of cracking in the vertical direction. 

Using an energy approach, it is possible to trace those components 
of the elemental stiffness matrix that are associated with each fracture code in 
every element in the LHR. This information can be assembled into the stiffness 
matrix such that a solution can be obtained for the incipient rupture of any 
subelement in the LHR. Note that this is most easily possible if (as in the 
present work), the constituent behavior is restricted to being linear elastic 
to the rupture point. Inelastic behavior, while more complicated, is not 
precluded, however. 


Interfacial 

Zone 




figure 3. DETAILED REPRESENTATION OF 
AN urn ELEMENT 




1 
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FIGURE 4. THE FRACTURE CODES CONTAINED 
IN AN LHR ELEMENT 
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Knowledge of the stiffness components attributable to each fracture 
condition also provides a method of calculating the energy-release rate for 
crack advance by any of the four codes provided for in each LHR element. Of 
course, the critical rupture energy values must be specified to provide a de- 
cision rule for breakage in each separate element. Note that while the model 
allows separate critical values to be specified for each of the four fracture 
codes for every one of the elements in the LHR, ordinarily different values 
will be specified only for the different constituents , i.e., fiber, matrix, or 
interface. 

Elemental Stiffness Formulation 


Energy principles can conveniently be used to determine the stiffness 
of the LHR structure. The sign conventions for an LHR element will be taken as 
shown in Figure 5. Then, using the following notation: 


til 

= rotational spring constant at the i — node 
K^^ = extensional stiffness in the x direction between the 


.th j .th , 

1 — and j — nodes 


K 




yy = extensional stiffness in the y direction between the 

.th j . th j 

X — and j — nodes 


= cross extensional component (effect of Poisson* s ratio) 
between the i— and j— nodes 

u^ = X displacement at the i^^ node 


v^ = y displacement at the i — node, 

the total strain energy U stored in the element for any set of arbitrary nodal 
displacements u^ and v^ can be written as 
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" ■ 1 <"2 - “!>' ^ 1 




2 XX 3 4 xy / 2 


3 “4) 




, 1 „23 , ^2 , „23 i^3 ^4 , ^2 '^l , , 

+ 2 Sy <''3-''2> + >'.<y ^ <''3-''2> 


, i 0 , 


2 1 Ax 


Ay 


+ 1 C ] ^^3 ^2^ 


2-2 Ay 


Ax 


1 ) ^‘"3“'"2^ ^V^^4^ 

+ ■# C„ ■( — T — — + 


2 3 Ay 


Ax 


+ 1 r , ^" 4 -^ 3 ) 

2 4 ) Ay Ax 


(7) 


It is desirable to determine an elemental stiffness matrix such that 


K 

e 



u 


A 
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This is accomplished by relating the node forces to the partial derivatives of 
the total strain energy- That is 


and 



F 


au 

9v. 

X 


(9) 


For example, 



( 10 ) 


Expansion of this relation will give the elements in the first row of the matrix 

[K ] . Subsequent derivatives taken with respect to the other degrees of freedom 
e 

generate the remaining rows. 

The spring constants used in the formulation must be related to the 
constituent’s elastic properties. This is accomplished via the relations 


K 


K 


■Ij ^ 

xy 

<^12 


ij = 

Qll 


XX 

Ax 

= 

n 

Ax 

yy 

CM 

CM 

r 

Ay 


" 4 *^33 


for i, j = 1, . . .4 
for i, j = 1, .... ,4 
for all i, j = 1, ... ,4 
for all i = 1, ... ,4 
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where [Q] is defined as the constitutive stiffness matrix in the usual manner 
as 


- - 


- - 


•- - 

a 




£ 

X 




X 

a 

= 

Q 


£ 

y 




y 

T 




Y 

xy 




xy 

- 


- 


„ _ 


( 11 ) 


Specifying these relations is tantamount to assuming that the material in each 

element is completely homogeneous. The desired heterogenity arises from the fact 

that the Q..’s are different for the fiber, matrix, and interface elements, 
ij 

The Energy-Release Rate 

The local rupture criterion used in this work is a generalized version 
of the ordinary energy-release rate (or crack-driving force) quantity of LEFM. 

In formulating it, attention must be placed on the complete structure. The 
discretized form of the equations of equilibrium for a complete structure can 
be represented by the matrix equation 


u = F , (12) 

where k is the structural stiffness matrix, u is the nodal displacement vector 
and F the applied nodal load vector. The energy-release rate G can be defined 
in a global sense by relating it to the change in work done by the applied load- 
ing and the strain energy during a virtual crack extension. This is 


k 


^ ~ da da 


( 13 ) 


where 


T 

W H u F 
^ 0 / 


( 14 ) 
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is the work done by the applied loading, 


u = i ( 15 ) 

is the strain energy and a represents the crack length. (In these equations, the 
superscript T denotes the transpose.) 

When Equations (14) and (15) are introduced into Equation (13), it 
is found that 


du 

"= di 




rr. dF 

T 

da 


1 T d ri n 

2 ^ da 


(16) 


It can be seen from Equation (12) that the quantity within the braces of Equa- 
tion (16) is a null vector. Furthermore, if is independent of the crack length, 
then the second term also vanishes. Equation (16) therefore reduces to 


G 


1 T d 

2 ^ da 


[k] u 
‘-a.-’ ^ 


(17) 


Now, it can be considered that the only contribution to the matrix ^ [k] is that 
due to the stiffness matrix K of the element containing the crack tip. If 
denotes the displacement vector of the nodes of the element only, then Equation 
(17) can be used to obtain an approximation of G that can be written as 


G = 


u*^ [1C ] u - u'^ ] u 

^ b ^ a ^ 

2Aa 


(18) 


where [K^] elemental stiffness matrix before (after) an extension 

Aa of the crack within the element. Recognizing that 


U 

e 


2 ^ 


[K] 


(19) 
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is the elemental strain energy, then an alternative form of Equation (18) is 


6 = 



( 20 ) 


A similar development for the energy-release rate occurs in the recently pub- 
lished papers by Hellen, et al. [6,7]« 


Program Solution Techniques 


To use the analysis that has been developed thus far, the material 
properties of the constituents and the micromechanical geometry must be specified. 
(The material properties are also used to determine the homogeneous orthotropic 
properties of the bulk composite through an effective modulus or volume fraction 
technique.) The way in which the constituent properties and the geometry are 
used to generate elemental stiffnesses for assemblage into the LHR will be des- 
cribed next. 

There are two classes of nodal points used in the LHR. The first 
consists of the points on the boundary of the LHR where displacements are pre- 
scribed. The nodal points in the interior of the LHR having unknown displacements 
make up the second class. It is useful to assemble the stiffness matrix in 
such a way that these two groups are isolated. The partition used here is 


’^22 

[ 1 

PI 

II 

r 

U) 



[ "33 1 


ll.i) 


In Equation (21), the subscripts do not refer to nodes, but to the nodal point 
classification just described. That is (u^) represents the prescribed displace- 
ment vector of the peripheral points while (u^) is the vector of the unknown dis- 
placements for which a solution is desired. Note that the number 1 was not used 
as a subscript because it is more convenient to reserve it for specified zero 
displacements . 

By performing the matrix multiplication indicated in Equation (21), 
the following two equations are obtained: 
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( 22 ) 


^32] I ^2 f [^33] 


{"3}={^3l= {0} 


( 23 ) 


Most often, it is either impossible or inefficient to store large matrices such 
as [^<^33] in the computer’s central processor core area. For this reason [1^33] 
is stored externally. The external storage can be minimized by utilizing the 
fact that the [^^3] matrix is always symmetric and banded; the latter term 
indicating the tendency for nonzero values to cluster around the diagonal of 
the matrix. By taking advantage of matrix partitioning (i.e., the grouping 
of prescribed and unprescribed degrees of freedom), symmetry, and the banded 
nature of the stiffness matrix, storage of only a small portion of the original 
total stiffness matrix is required. However, this reduced storage requirement 
can still be very large in some problems. 

Solution of Equation ( 23 ) is accomplished through the use of 
Gaussian elimination with back substitution. Small parts of stored [^33] matrix 
are brought into core as required in this operation. Assemblage of the 

■k 

required parts of the stiffness matrix is done by the direct stiffness method. 
Components that would ordinarily fall within the [^ 22 ^ ^^ 23 ^ submatrices 

are discarded. Stiffness components that fall within the [k32l submatrix are 
used in a nonzero vector. Finally, those values of [^33] falling on or 

above its diagonal are stored externally. 

Once the boundary conditions have been prescribed, an out-of-core 
type solution of the partitioned LHR stiffness matrix with boundary conditions 
for a unit load is performed. The resulting displacement vector is then used 
to calculate the potential fracture energies associated with each of the four 
fracture codes. Ratios of the prescribed critical energy levels to these cal- 
culated energies are next computed. The applied load is then adjusted such 
that the highest ratio becomes exactly equal to unity. This critical region 
(if one exists) is allowed to break in one of the four ways (codes) described 
above. Appropriate modifications are then made to the LHR stiffness matrix. 

* A detailed description of the direct stiffness assemblage technique used 
in this work is given by Desai and Abel [8]. 
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Because the system is completely linear elastic, a solution for another load level 
can then be performed in exactly the same manner. In this way, the properties 
of a crack tip damage zone as a function of an increasing applied load can be 
generated. 

As each local rupture occurs, the prescribed boundary conditions must 
be adjusted- If the crack has propagated in a self-similar fashion, this is 
readily done by shifting the origin (cf, Figure 2). If not, this may present 
some very real difficulties for then the continuum solution no longer exactly 
applies and some approximations must be used for setting the boundary conditions. 
The larger the LHR, the less likely this will be necessary. Of course, large 
LHR’s require large solution times even though storage is generally not a problem 
because of the out-of-core solution. It is likely that the boundary condition 
problem can be handled by introducing a flexible boundary-condition scheme (see 
Recommended Future Research section), or, possibly, by allowing the fractured LHR 
to interact in a hybrid continuum model for future boundary updating. 

Verification of Computational Mode l 

Before performing extensive calculations on composite materials, it 
is appropriate to verify the computational model by checking it against known 
solutions. By letting the model represent a linear elastic homogeneous material, 
two checks can be made. These are on the displacements in the LHR and on the 
calculated energy-release rate. 

Consider the element configuration shown in Figure 1- A check on 
the calculated displacements can be made by imposing the continuum derived 
boundary conditions on a LHR region where the LHR is used to simulate a com- 
pletely homogeneous region. Displacements at the node points were calculated 
and used to calculate average stresses within the elements. The computed 
values of stress and displacement were then compared with values calculated by 
the continuum solution for the interior of the region. Agreement between the 
two solutions was found to be quite good. The check thus provided one necessary 
verification of both the LHR element itself and of the finite element assemblage 
and solution procedure as well. 

As described in the above, all of the energy released within the 
element at the crack tip, if an incremental amount of crack extension were to 
occur, should equal the strain energy-release rate, G. To check this, cal- 
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culations were made using the exact near crack tip displacement field for an 
elastic isotropic material. The energy-release rate obtained from this calcula- 
tion is denoted by G * to distinguish it from the exact theoretical value G. The 
necessity for introducing this distinction is as follows. 

The quantity G can be formally defined in terms of a virtual crack 
extension as the total change in energy of the body per unit area of crack 
extension; cf. Equation (13). The quantity G*, on the other hand, reflects 
only the change in energy of the element near the crack tip. Consequently, G* 
does not include the change in energy of the remainder of the body as the crack 
extends. In order to account for the contribution arising from the change in 
energy of the remainder of the structure, a term denoted by 6 G * was formulated. 
The sum G* + 6G* is then taken as the appropriate approximation to G in this 
work. 

A numerical experiment was performed to calculate G* and 6G* as a 
function of the aspect ratio Ay/ Ax of the LHR elements. The results are shown 
in Table II. It is quite evident that the sum G* + 6G* provides an accurate 
approximation to G for the elongated aspect ratios that are convenient to use 
in the calculations on fiber composite materials. 

There are two further points of interest about the results shown in 

Table I. First, while it may be surprising that 6G* and G* are functions only 
of the ratio of Ay to Ax, this is a direct consequence of the use of rigid 
boundary conditions which put the near tip displacements into the energy rela- 
tionship used to calculate G. The resulting expression contains only the 
aspect ratio, Ay/Ax. Second, while it might be expected that 6G* should con- 
verge to zero as the grid spacing is collapsed (i.e., as both Ay ^ 0 and Ax->-0), 
it can be seen that this was not the case. Instead, convergence is obtained 
only as the ratio of Ay to Ax becomes large. Then, 6G* approaches a value equal 
to 11 percent of G while G* converges to value equal to 86 percent of G. 

A number of different computations for heterogeneous materials have 
also been made. These runs have been made with a model similar to the one shown 
in Figure 3. Preliminary results indicate that most of the initially envisioned 
modes of failure are exhibited by the model. A discussion of these results is 
given in the next section of this report. 
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TABLE I* AN EXAMPLE CALCULATION OF THE STRAIN ENERGY 

RELEASE RATE AS A FUNCTION OF THE GRID SIZE RATIO 
FOR A HOMOGENEOUS LINEAR ELASTIC MATERIAL 


M 

Ax 

G* 

G 

G 

(G* + dG*) 
G 

1 

0.71 

0.75 

1.46 

2 

0.72 

0.42 

1.14 

4 

0.77 

0.25 

1.02 

8 

0.81 

0.18 

0.99 

16 

0.83 

0.15 

0.98 

32 

0.84 

0.13 

0.97 

64 

0.85 

0.12 

0.97 

128 

0.86 

0.11 

0.97 

256 

0o86 

0,11 

0.97 
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EXAMPLE COMPUTATIONAL RESULTS FOR 
FRACTURE OF FIBER-REINFORCED COMPOSITES 


A number of computations have been performed using the analysis 

technique described in the preceding section of this report. As discussed 

in this section, the results demonstrate the ability of the model to exhibit 

most of the actual failure mechanisms in fiber-reinforced composite materials. 

These include fiber-matrix debonding, fiber bridging, matrix bridging, and 

matrix crazing. The only important micromechanical mechanism that the model 

* 

Currently does not explicitly represent is fiber pull-out. 

Calculations With Arbitrarily Varied Rupture Properties 

The following describes a series of example calculations in which 
the various fracture mechanisms are made to manifest themselves. The re- 
sults shown here were obtained by varying each of the constituent’s rupture 
properties and, in some instances, the constituent elastic properties, over a 
wide range of values. Hence, the actual values selected to perform these 
example computations are not totally realistic. The intent here is to provide 
a qualitative verification that the model is capable of achieving its primary 
function rather than to make realistic predictions. 

Unless otherwise stated, the elastic constants used in the calcula- 
tions are those given in Table ii. 

TABLE III. ELASTIC PROPERTIES USED IN THE SIMULATION 
OF A GRAPHITE EPOXY COMPOSITE 


Constituent 

E 

Elastic 
Modulus (ksi) 

V 

Poisson ' s 
Ratio 

Fiber 

28,000 

0.3 

Matrix 

495 

0.3 

Interface 

495 

0.3 


* As discussed in the Recommended Future Research section of this report, 
steps to incorporate both fiber pull-out and interply delamination can 
be performed in a straightforward manner in subsequent work. 
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These properties were intended to be nominally equivalent to those of a 
graphite epoxy composite with a volume fraction of fiber equal to 70 percent. 
Consistent with these values are the following properties of the bulk composites: 


= 18,000 ksi, 0.25 

= 690 ksi, G = 3,000,000. 


Note that the constituents were assumed to be isotropic. However, the model 
does not require this. 

The first computation to be discussed is shox%na in Figure 6. This 
result illustrates the fiber-matrix debonding mechanism, a mechanism that occurs 
in a large number of cases. Typically, a crack will advance through the matrix 
and then vertically separate or axially split the interface. In this example, 
the interfacial elements were made relatively weak while the fiber and matrix 
material were given relatively average strengths. The loading in this example 

;'c 

was purely Mode I. It was found that axial splitting is even more prevalent 
when Mode II loads are applied. 

Matrix crazing typically occurs when a crack advancing through the 
matrix reaches a very strong stiff fiber. Figures 7 and 8 show typical examples. 
If the fiber does not break, a number of local events occur in the matrix that 
seem to approximate matrix crazing. Typically, the crazing does not occur on 
the uncracked side of the fiber. The reason is that the displacements on the 
cracked side of the fiber are "locked in" by the adjacent highly stretched 
fiber. Consequently, crazing usually occurs when strong very stiff fiber 
properties are inserted into the model together with average strength and stiff- 
ness properties for the matrix and interface. 

Matrix bridging is shown in Figure 9. This phenomenon typically 
occurs when the fibers are considered to be stiff, but weak. In these in- 
stances the matrix and interface are capable of withstanding higher elonga- 
tion than the fiber. So, they will remain intact while the fiber cracks. In 
the example shown in Figure 9, some additional interfacial response can also 
be noticed. 


* Unless otherwise stated, in the calculations described in this section of the 
report, the applied loading was a tensile normal stress in the direction 
parallel to the fibers. 




FIGURE 6. EXAMPLE CALCULATION WITH WEAK 


EVENT # 

COMPONENT 

RELATIVE LOAD 

1 

MA.TRIX 

1.0 

2 

MATRIX 

1.0 

3 

INTERFACE 

2.2 

4 

INTERFACE 

4.9 

5 

INTERFACE 

4.9 

6 

INTERFACE 

5.3 

7 

INTERFACE 

5.3 

8 

INTERFACE 

7.5 

9 

INTERFACE 

7.5 

10 

INTERFACE 

7.8 

11 

INTERFACE 

7.8 
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M 
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FIGURE 7. EXAMPLE CALCULATION WITH STRONG STIFF FIBERS AND WEAK INTERFACE SHOWING MATRIX CRAZING 






EVENT # 



COMPONENT 

RELATIVE LOAD 

muRix 

1.0 

MATRIX 

6.1 

INTERFACE 

7.3 

MATRIX 

27.5 

MATRIX 

27.5 

MATRIX 

27.5 

INTERFACE 

27.5 

MATRIX 

27.8 

MATRIX 

27.8 

MATRIX 

27.8 

INTERFACE 

27.8 

MATRIX 

39.5 











FIGURE 9, EXAMPLE CALCULATION WITH STIFF WEAK FIBERS SHOWING MATRIX BRIDGING 
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Fiber bridging is more prevalent when the fiber modulus does not 
greatly exceed the matrix modulus. (In graphite epoxy, the modulus ratio is 
on the order of 50 and fiber bridging does not then seem to occur.) Figure 10 
demonstrates this behavior. Fiber bridging can also be related to the micro- 
structural geometry. That is, in computations performed with relatively thin 
fibers, this mechanism becomes more dominate. 

It should be emphasized that the examples given in the above deliber- 
ately used exactly the same LHR configuration and applied loading. The program 
has the flexibility to consider different loading conditions and LHR geometries. 
However, this option was not exercised here because of the complications that 
would be added to the interpretation of the results. That is, it would be 
difficult to separate the effects due to the LHR geometry and loading from 
those due to variations in the rupture strength and modulus. Such calculations 
will be deferred until further refinements have been incorporated into the 
model. Again, the purpose of the computations given in Figures 6-10 is to give 
a qualitative demonstration that the model is capable of coping with the micro- 
mechanical failure processes involved in the fracture of composite materials, 
not to produce precise quantitative results. 

Finally, note that the relative load levels at which each individual 
fracture event occurs is recorded. (These levels are relative to the load 
level at which the initial rupture event occurs.) It can be seen that the 
load level must ordinarily be increased to obtain additional rupture events, 
or, in other words, in order to enlarge and propagate the crack-tip damage 
zone. This can be contrasted with the behavior of completely homogeneous, 
perfectly brittle materials represented by LEFM where crack extension, once 
initiated, would continue iri a catastrophic manner under a constant load level. 
The initial stages of the fracture event in fiber composite materials therefore 
can obviously be characterized in terms of a stable growth process. This sug- 
gests that a modification of the crack growth resistance curve (R curve) ap- 
proach developed for the fracture of ductile materials could be useful for 
studying fractures in fiber composite materials. Such an approach, it might 
be mentioned, has already been suggested in the literature; for example, by 
Gaggar and Broutman [9], However, these approaches are usually semiempirical 
in nature, relying heavily on LEFM concepts. The present approach, in con- 
trast, should be able to make a direct prediction of the crack growth resis- 
tance parameter purely from fundamental-level considerations. 
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EVENT 

COMPONENT 

RELATIVE LOAD 

1 

MATRIX 

1. 

2 

MATRIX 

1.0 

3 

INTERFACE 

1.6 

4 

INTERFACE 

1.6 

5 

MATRIX 

7.5 

6 

MATRIX 

7.5 

7 

INTERFACE 

7.5 

8 

INTERFACE 

7.5 

9 

INTERFACE 

7.5 

10 

INTERFACE 

7.5 


STIFF MATRIX SHOWING FIBER BRIDGING 
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Comparison of Calculated Results 
with VPISU Experimental Data 

The computations described in the preceding section demonstrate 
that the model developed in this report displays the various micromechanical 
failure mechanisms actually exhibited by fiber reinforced composites. A 
more quantitative verification is also possible. This can be done by com- 
paring the predictions of the model with the experimental results obtained 
in the concurrent NASA-Ames sponsored research at Virginia Polytechnic 
Institute and State University (VPISU) by Brinson and Yeow [10]. They have 
measured the strengths of both unidirectional and angle ply graphite/epoxy 
composites using unnotched, single edge notched, and double edge notched speci- 
mens with the crack introduced at various angles to the fiber direction. In 
this section, the model will be used to estimate the strengths of some of the 
Brinson-Yeow unidirectional notched tests using input values inferred from 
their unnotched tests. 

A precise prediction of failure loads is not to be expected at the 
present stage of this research for several reasons. First, because of the 
stable damage growth that precedes fracture in composites, computations per- 
formed using rigid boundary conditions (which artifically constrain the damage 
growth process) will not be realistic. The predictions of the initial local 
failure event — the threshold of stable damage growth — should be reasonably 
well predicted, however. This load level will therefore be taken as the 
prediction of the model to be compared with the experimental measurements . 

It should be recognized that, because of the neglect of the stable growth 
regime, these should underestimate the actual failure loads. 

A second reason for the lack of precision inherent in the present 
predictive capability is that the constituent properties needed in the model 
have not been properly determined as yet. As in the preceding section, linear 
elastic-perf ectly brittle behavior can be assumed with handbook values being 
taken for the elastic properties. Rupture properties are not as readily 
available, however. To circumvent this difficulty, the experimental results 
on unnotched specimens obtained by Brinson and Yeow can be used. Their data 
on the strengths of unnotched coupons pulled to failure with the load in the 
fiber direction (0 = 0°) and normal to the fiber direction (0 = 90°) are given 
in Table III. Relying on their observation that matrix failure occurred in 
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TABLE III. CRITICAL VALUES FOR MATRIX FAILURE 
IN GRAPHITE EPOXY COMPOSITES 


Angle 

Between 


(a) 

Experimental Values 


Critical 

Load and 

Elastic 

Fracture 

Fracture 

Strain Energy 

Fiber 

Modulus 

Stress 

Strain 

Density 

Directions 

(ksi) 

(ksi) 

(%) 

(in. lb/ in. 3) 

0° 

18,200 

154.7 

0.82 

576 

o 

O 

1,700 

6.1 

0.44 

13.6 


(a) Data of Brinson and Yeow [10] on unidirectional unnotched specimens for 
a loading rate of 2 x 10"^ inches/min. 
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virtually all cases, the critical strain energy density quantity can be cal- 
1 2 

culated (i.e., ^ model. These figures are given in the 

last column in Table III. 

As discussed on page 11 of this report, the applied loads acting 
on the body are communicated to the crack tip elements via the anisotropic 
elastic stress intensity factors. F.r test specimens, such as the single and 
double edge notched configurations used by Brinson and Yeow, there is a term 
in the stress intensity factor that contains the effects of the finite geometry. 
This term is not the same for anisotropic and isot . opic bodies. However, in 
view of the approximate nature of the present calculations, the refinement 
added through the use of the rather complicated anisotropic expression was 
not believed to be warranted. Hence, the isotropic expressions were used. 

For Mode I loading, these have the form. 


1 

= a Y (24) 

where a is the applied stress normal to the crack plane, a is the crack length 
and Y is a dimensionless function of the specimen geometry. For the double 
edge notched configuration 


Y = 1.99 + 0.76 (-^) - 8.48 (-^)^ + 27.36 , (25) 

while for the single edge notched configuration 

Y = 1.99-0.4l(^) + 18.7 - 38.48 + 53.85 (26) 

where W is the plate width. 

For a given crack length and load-fiber orientation in either the 
single or double edge notched configuration, the load corresponding to the 
threshold of damage is determined by a single calculation. This is due to 
the linear elastic material behavior assumed in the current model. That is, 
a computation can be performed for an estimated value of K^. The solution 
can then be scanned to determine the highest strain energy density in any 
matrix element. By "scaling up" the solution to force this value to match 
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the critical value given in Table IV, the value of corresponding to the 
threshold of damage is determined. The final step is to use Equation (24) 
to calculate the applied stress for direct comparison with the experimental 
results . 

Comparisons of the calculated results as a function of crack length 
with the Brinson-Yeow results on unidirectional composites are shown in 
Figures 11 and 12 for the single and double edge notch configurations, re- 
spectively. Note that a semilog format is used (so that both the 0° and the 
90° fiber-load angle results can be shown on the same plot) and this may make 
the agreement appear to be better than it really is. Nevertheless, in view of 
the remarks made in the preceding section, the prediction is reasonably 
accurate and, as expected, generally provides a lower bound to the experimental 
data. Note also that the unnotched result (-^ = 0) is included on the theo- 
retical curve to emphasize that this point was used to construct the curve. 

Finally, comparisons with experimental data can also be made for the 
case of cracks introduced at an angle to the fiber direction. An example 
computation is shown in Figure 13. The threshold-of-damage load level calcu- 
lations are made as above except that both the Mode I and Mode II stress 
intensity factors now are involved in the computation. These have been taken 
as 


- 2 
Kj = aa^ Y sin 

and 

1 

= aa^' Y cos <i> sin c() , (27) 

where ()> is the angle between the crack plane and the fiber direction. A 
comparison between the predicted results and the Brinson-Yeow experimental 
results as a function of <}) for a 0° load-fiber angle is shown in Figure 14. 
Again, it can be seen that the model provides a reasonable lower bound to 
the experimental results. 



Remote Applied Stress, ksi 





Remote Applied Stress , ksi 










Remote Applied Stress . ksi 



denote failure loads determined at 
VPISU on unidirectional graphite 
epoxy composites for zero angle 
between, the load direction and 
the fiber direction 


Angle Between Crack Plane and 
Load Direction (Degrees) 


FIGURE 14, COMPARISON OF PREDICTED STABLE GROWTH THRESHOLD 
WITH EXPERIMENTAL FAILURE LOADS FOR SINGLE EDGE 
NOTCH SPECIMENS WITH ANGLE CRACKS 
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RECOMMENDED FURTHER RESEARCH 


Being two dimensional, the model described in this report is so far 
limited in application to unidirectional laminates and can only reflect the 
local damage modes that occur in plane deformation. In subsequent work, it 
would be desirable to extend the model from two dimensions to three dimensions 
to treat angle ply laminates. In addition to treating the effect of different 
layups, new modes of damage (e,g., interply delamination, fiber pull-out) 
should also be included. 

A straighforward extension of the model could be made by devising 
a three-dimensional LHR. But, the number of storage locations required for 
the three-dimensional grid would likely exceed the fast access storage of most 
computers. To circumvent this situation, the present two-dimensional model 
could be extended to angle ply laminates by considering laminates with an LHR 
only in a single ply. Figure 15 shows this concept. In the simplest case, 
the LHR can be confined to the most critical ply as the load level is increased. 
In general, there can be an LHR in each ply with attention shifted from one ply 
to another in turn. In either case, it would be necessary to assume that the 
interply interactions are not affected by the precise details of the local 
damage occurring in each ply. 

In Figure 15a, a laminate containing a flaw is shown with generalized 
loading conditions. In this the first level of modeling, all of the plys are 
considered as homogeneous orthotropic layers. Displacements can be calculated 
in this model for any load level and applied to the boundaries of the local 
homogeneous region in each ply. Then, just as in the two-dimensional model, 
the LHR shown in Figure 15b can be monitored to determine the local rupture 
events and the extent of the damage zone at that load level. When significant 
crack growth has occurred in one or more plys, the first level model must, of 
course, be made to reflect this. But, the basic approach would not be changed 
thereby. 

Two distinct techniques must be developed in order to obtain an 
accurate set of boundary displacements for each LHR when extensive damage 
occurs in that or in neighboring plys. The first is to incorporate the effect 
of the stiffness change in the LHR into the first level model of the laminate. 
The loading conditions would then be obtained from this modified first level 
model (i.e., one with locally reduced stiffnesses) to give the new boundary 
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conditions for the LHR. In the second technique, the influence of the stiff- 
ness of the laminate away from the crack tip will be directly transferred to 
the LHR through what can be called a "flexible boundary condition" technique. 

This approach can be likened to placing a series of springs on the LHR boundary 
with stiffnesses taken from the first level model. The actual technique is 
more sophisticated than this, however, as follows. 

In the work described so far, the peripheral elements in the LHR 
have been positioned according to the displacement field given by a completely 
linear elastic continuum solution for the given crack tip location and applied 
loads. This is not strictly correct even when the displacements are periodically 
updated as the crack extends in the LHR. Even in the absence of local damage, 
the highly nonlinear inhomogeneous nature of the crack-tip region in a com- 
posite will cause significant departures from the continuum displacements. A 
scheme for estimating these displacements based on a determination of the equi- 
librium state between the LHR and the continuum can be obtained by adapting 
the technique reported in Reference 3. Such a scheme has been referred 
to as the "flexible boundary condition" approach. The approach used so far can 
be called the "rigid boundary condition" approach because, even with periodic 
updating to reflect the progress of the crack through the LHR, the LHR does 
not directly affect the continuum region surrounding it. The flexible boundary 
condition approach, in contrast, accounts for the interaction. 

Finally, the range of inodes of damage should be extended over those 
of the two-dimensional model to include interply delamination and a "free- 
edge" effect. The general approach could be similar to that described in 
References 10 and 11. Attention should first be focused on through-the-thickness 
cracks in laminated plates under tension with part-through flaws being considered 
later. Of most importance, the three-dimensional model of an angle ply laminate 
should permit arbitrarily varied multidirectional layups to be explicitly con- 
sidered. The properties for the ply moduli could be obtained from laboratory 
investigations, e.g., in the program being conducted for NASA-Ames at Virginia 
Polytechnic Institute and State University. 
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